#install.packages("tidyverse")
#install.packages("readstata13")
#install.packages("interflex")

rm(list=ls(all=TRUE))

library(tidyverse)
library(readstata13)
library(interflex)

# Set Working Directory
#setwd("")

dat<-read.dta13("FDI&Conflict.dta")

################################################################################
###### Figure A: Plot of Untransformed and Transformed FDI Variables ###########
################################################################################

pdf(file="fdi_trans.pdf",width=6.5,height=2.2)
par(mfrow=c(1,3),mar=c(2.1,2.2,.5,.8),xpd=F,
    cex.lab=.6, cex.axis=.6,cex=.7,mgp = c(1.15,.5,0))
hist(dat$rfdicap/1000,breaks=50,col="grey80",main="",
     xlab="Real FDI per Capita in Thousands of $")
hist(dat$cr_rfdicap,breaks=50,col="grey80",main="",
     xlab="Real FDI per Capita (Cube Root)")
hist(log(dat$rfdicap+8292),breaks=50,col="grey80",main="",
     xlab="Real FDI per Capita (log)")
dev.off()


################################################################################
######################### Model 1 in Table 4: RPE ##############################
################################################################################
mydata1<-dat%>%select(ccode, year, cr_rfdicap, onset2v414, incidencev414, lwdist,
                     lmrpe, llrgdpcap, llpop, llgrowthrate, lpolity2, 
                     lloil_gas_cap, elf85, relfrac, lmtnest, ncontig, coldwar, 
                     t1_onset2, t2_onset2, t3_onset2)%>%drop_na()

mydata1<-mydata1%>%mutate(s = as.numeric(incidencev414==1&onset2v414==0))%>%
                filter(s==0)

fit.fdi1<-lm(cr_rfdicap~lwdist+lmrpe+llrgdpcap+llpop+llgrowthrate +lpolity2 +
               lloil_gas_cap+ elf85+relfrac+ lmtnest+ ncontig+ coldwar+ 
               t1_onset2+ t2_onset2+ t3_onset2, data=mydata1)

mydata1$fdi.hat1<-predict(fit.fdi1)

fit.fdi1a<-lm(fdi.hat1~lmrpe+llrgdpcap+llpop+llgrowthrate +lpolity2 +
                lloil_gas_cap+ elf85+relfrac+ lmtnest+ ncontig+ coldwar+ 
                t1_onset2+ t2_onset2+ t3_onset2, data=mydata1)
mydata1$fdi.res1<-resid(fit.fdi1a)

fit.onset1<-lm(onset2v414~lmrpe+llrgdpcap+llpop+llgrowthrate +lpolity2 +
                 lloil_gas_cap + elf85+relfrac+ lmtnest+ ncontig+ coldwar+ 
                 t1_onset2+ t2_onset2+ t3_onset2, data=mydata1)

mydata1$onset.res1<-resid(fit.onset1)

# Common Support
out.raw.rpe<-interflex(estimator = "raw", Y = "onset.res1", D = "fdi.res1", 
            X = "lmrpe", data = mydata1, Ylabel="Probability of Conflict Onset", 
            Dlabel="Real FDI per Capita (Predicted)", Xlabel="RPE",ncols=3)

pdf(file="commonsupport_rpe.pdf",
    width=6.25,height=3)
plot(out.raw.rpe)
dev.off()

# Binning Plot
covars<-c("llrgdpcap","llpop","llgrowthrate","lpolity2","lloil_gas_cap","elf85",
   "relfrac","lmtnest","ncontig","coldwar","t1_onset2","t2_onset2","t3_onset2")

out.binning.rpe <- interflex(estimator = "binning", Y = "onset2v414", D ="fdi.hat1", 
                X ="lmrpe", Z =covars, data=mydata1, vcov.type  = "robust", 
                Xlabel = "Relative Political Extraction", Dlabel = "FDI", 
                Ylabel = "Conflict Onset", cex.lab = 1)

pdf(file="binning_rpe.pdf",width=6.25,height=5)
plot(out.binning.rpe)
dev.off()

# Wald Test
print(out.binning.rpe$tests$p.wald)

#Kernel Regression
set.seed(123)
out.kernel.rpe<-interflex(estimator = "kernel", Y = "onset2v414", D = "fdi.hat1",
                X = "lmrpe", Z = covars, data = mydata1, nboots = 200, 
                parallel = TRUE,  Ylabel="Conflict Onset", Dlabel="FDI", 
                Xlabel="Relative Political Extraction")

pdf(file="kernel_rpe.pdf",width=6.25,height=5)
plot(out.kernel.rpe)
dev.off()

################################################################################
####################### Model 2 in Table 4: Tax/GDP ############################
################################################################################
mydata2<-dat%>%select(ccode, year, cr_rfdicap, onset2v414, incidencev414, lwdist,
                      lltax_m, llrgdpcap, llpop, llgrowthrate, lpolity2, 
                      lloil_gas_cap, elf85, relfrac, lmtnest, ncontig, coldwar, 
                      t1_onset2, t2_onset2, t3_onset2)%>%drop_na()

mydata2<-mydata2%>%mutate(s = as.numeric(incidencev414==1&onset2v414==0))%>%
  filter(s==0)

fit.fdi2<-lm(cr_rfdicap~lwdist+lltax_m+llrgdpcap+llpop+llgrowthrate +lpolity2 +
               lloil_gas_cap+ elf85+relfrac+ lmtnest+ ncontig+ coldwar+ 
               t1_onset2 + t2_onset2 + t3_onset2, data=mydata2)

mydata2$fdi.hat2<-predict(fit.fdi2)

fit.fdi2a<-lm(fdi.hat2~lltax_m+llrgdpcap+llpop+llgrowthrate +lpolity2 +
                lloil_gas_cap+ elf85+relfrac+ lmtnest+ ncontig+ coldwar+ 
                t1_onset2 + t2_onset2 + t3_onset2, data=mydata2)

mydata2$fdi.res2<-resid(fit.fdi2a)

fit.onset2<-lm(onset2v414~lltax_m+llrgdpcap+llpop+llgrowthrate +lpolity2 +
                 lloil_gas_cap + elf85+relfrac+ lmtnest+ ncontig+ coldwar+ 
                 t1_onset2 + t2_onset2 + t3_onset2, data=mydata2)

mydata2$onset.res2<-resid(fit.onset2)

# Common Support
out.raw.tax<-interflex(estimator = "raw",Y = "onset.res2", D = "fdi.res2", 
             X = "lltax_m", data = mydata2, Ylabel="Conflict Onset", 
          Dlabel="Real FDI per Capita (Predicted)", Xlabel="Tax/GDP",ncols=3)

pdf(file="commonsupport_tax.pdf", width=6.25,height=3)
plot(out.raw.tax)
dev.off()

# Binning Plot
out.binning.tax <-interflex(estimator = "binning", Y = "onset2v414", 
                  D ="fdi.hat2", X ="lltax_m", Z =covars, data=mydata2, 
                  vcov.type  = "robust", Xlabel = "Relative Political Extraction",
                  Dlabel = "FDI", Ylabel = "Conflict Onset", cex.lab = 1)

pdf(file="binning_tax.pdf",width=6.25,height=5)
plot(out.binning.tax)
dev.off()

# Wald Test
print(out.binning.tax$tests$p.wald)

# Kernel Regression
set.seed(123)
out.kernel.tax <- interflex(estimator = "kernel", Y = "onset2v414", 
                  D = "fdi.hat2", X = "lltax_m", Z = covars, data = mydata2, 
                  nboots = 200, parallel = TRUE, Ylabel="Conflict Onset",  
                  Dlabel="FDI", Xlabel="Tax (% of GDP, log)" )

pdf(file="kernel_tax.pdf",width=6.25,height=5)
plot(out.kernel.tax)
dev.off()

################################################################################
############# Model 3 in Table 4: Primary School Enrollment Rate ###############
################################################################################
mydata3<-dat%>%select(ccode, year, cr_rfdicap, onset2v414, incidencev414, lwdist,
                      lenrollpri_m, llrgdpcap, llpop, llgrowthrate, lpolity2, 
                      lloil_gas_cap, elf85, relfrac, lmtnest, ncontig, coldwar, 
                      t1_onset2, t2_onset2, t3_onset2)%>%drop_na()

mydata3<-mydata3%>%mutate(s = as.numeric(incidencev414==1&onset2v414==0))%>%
  filter(s==0)

fit.fdi3<-lm(cr_rfdicap~lwdist+lenrollpri_m+llrgdpcap+llpop+llgrowthrate+
            lpolity2 +lloil_gas_cap + elf85+relfrac+ lmtnest+ ncontig+ coldwar+ 
            t1_onset2 + t2_onset2 + t3_onset2, data=mydata3)

mydata3$fdi.hat3<-predict(fit.fdi3)

fit.fdi3a<-lm(fdi.hat3~lenrollpri_m+llrgdpcap+llpop+llgrowthrate +lpolity2 +
                lloil_gas_cap+ elf85+relfrac+ lmtnest+ ncontig+ coldwar+ 
                t1_onset2 + t2_onset2 + t3_onset2, data=mydata3)
mydata3$fdi.res3<-resid(fit.fdi3a)

fit.onset3<-lm(onset2v414~lenrollpri_m+llrgdpcap+llpop+llgrowthrate +lpolity2 +
                 lloil_gas_cap+ elf85+relfrac+ lmtnest+ ncontig+ coldwar+ 
                 t1_onset2 + t2_onset2 + t3_onset2, data=mydata3)

mydata3$onset.res3<-resid(fit.onset3)

# Common Support
out.raw.enroll<- interflex(estimator = "raw", Y = "onset.res3", D = "fdi.res3", 
                 X = "lenrollpri_m", data = mydata3, 
                 Ylabel="Probability of Conflict Onset", 
                 Dlabel="Real FDI per Capita (Predicted)", Xlabel="Enrollment", 
                 ncols=3)

pdf(file="commonsupport_enroll.pdf",width=6.25,height=3)
plot(out.raw.enroll)
dev.off()

# Binning Plot
out.binning.enroll <- interflex(estimator = "binning", Y = "onset2v414", 
                      D ="fdi.hat3", X ="lenrollpri_m", Z =covars, data=mydata3, 
                      vcov.type  = "robust", Xlabel = "Primary School Enrollment Rate",
                      Dlabel = "FDI", Ylabel = "Conflict Onset", cex.lab = 1)

pdf(file="binning_enroll.pdf",width=6.25,height=5)
plot(out.binning.enroll)
dev.off()

# Wald Test
print(out.binning.enroll$tests$p.wald)

#Kernel Regression
set.seed(123)
out.kernel.enroll<-interflex(estimator = "kernel", Y = "onset2v414", D = "fdi.hat3", 
                   X = "lenrollpri_m", Z = covars, data = mydata3, nboots = 200, 
                   parallel = TRUE, Ylabel="Conflict Onset", Dlabel="FDI", 
                   Xlabel="Primary School Enrollment Rate" )

pdf(file="kernel_enroll.pdf",width=6.25,height=5)
plot(out.kernel.enroll)
dev.off()
